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An Asymptotic Preserving Two-Dimensional Staggered Grid Method 

for multiscale transport equations 
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Abstract 

We propose a two-dimensional asymptotic preserving scheme for linear transport equations with 
diffusive scalings. It is an extension of the time splitting developed by Jin, Pareschi and Toscani [15] , 
but uses spatial discretizations on staggered grids, which preserves the discrete diffusion limit 
with a more compact stencil. The first novelty of this paper is that we propose a staggering in 
two dimensions that requires fewer unknowns than one could have naively expected. The second 
contribution of this paper is that we rigorously analyze the scheme of Jin, Pareschi, and Toscani 
m- We show that the scheme is AP and obtain an explicit CFL condition, which couples a 
hyperbolic and a parabolic condition. This type of condition is common for asymptotic preserving 
schemes and guarantees uniform stability with respect to the mean free path. In addition, we 
obtain an upper bound on the relaxation parameter, which is the crucial parameter of the used 
time discretization. Several numerical examples are provided to verify the accuracy and asymptotic 
property of the scheme. 


1 Introduction 


The linear transport equation models particles interacting with a background medium (e.g., neutron 
transport, linear radiative transfer). In general, the model in scaled variables can be written as [19] 


ed t f + v • V x / 


Oj_ 
2 7T 




+ sQ, 


(1) 


where /(i,x, u) denotes the probability density distribution depending on time t, position x £ M 2 , 
and direction of velocity v = (£, rj) & il = {(£,??) : —1 < £,r? < 1, £ 2 + r/ 2 = 1} (two-dimensional 
flatland model - the extension to three dimensions is straightforward). Moreover, oy = a s +e 2 cr a is the 
total transport coefficient, a s is the scattering coefficient, cr a is the absorption coefficient, and Q is a 
u-independent source term. It is well-known that the limiting equation (e —> 0) of 0 is the diffusion 
equation: 

d t p = jVx ' - a a p + Q, (2) 

where p(t, x) = A fn f{t, x, v)dv. 
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In many applications, the scaling parameter of the transport equation e (mean free path) may differ 
in several orders of magnitude, ranging from the rarefied kinetic regime to the hydrodynamic diffusive 
regime. When e is small, in the diffusive regime, the equation becomes numerically stiff, which leads to 
numerical challenges: straightforward explicit implementations lead to high computational costs in the 
diffusive regimes; fully implicit schemes could be difficult to implement [6j; multiscale multiphysics 
domain-decomposition approaches, which couple models at different scales, have difficulties in the 
transition zones, since they need to transfer data from one scale to another P33- Thus, it is desirable 
to develop schemes which are suitable for all regimes (no domain-decomposition), but do not require 
a resolved grid in space and small time compared to the mean free path. This is the objective of 
asymptotic preserving (AP) schemes. 

A scheme is called AP if it preserves the discrete analogue of the asymptotic transition from the 
microscopic scale to the mascroscopic one 13 IE], namely, in the limit e —> 0, the discretization of the 
above transport equation ([!]) should yield a discretization of the diffusion equation ([2]). Such schemes 
allow mesh sizes and time steps much bigger than the mean free path/time, yet still capture the 
correct physical behavior. The development of such schemes started with steady problems of linear 
transport equations by Larsen, Morel, and Miller [23] and for boundary value problems by Jin and 
Levermore mm- Uniform convergence with respect to the mean free path for an AP scheme was 
first established by Golse, Jin, and Levermore [9] . 

In [23], and its follow-up [22] . several space discretizations for steady transport problems were 
investigated, among them diamond differencing and a linear discontinuous Galerkin (LD) method, 
both of which were identified to behave well in the asymptotic regime. Furthermore, LD leads to a 
compact stencil for the diffusion equation. 

For time-dependent problems, AP schemes were first designed for nonlinear hyperbolic systems 
with relaxation by Jin and Levermore m IT5l . There one needs to design both the time and the 
spatial discretization carefully m , in particular, to overcome the stiffness of the source term. AP 
schemes for time-dependent transport equations with diffusive scaling started by Jin, Pareschi, and 
Toscani nans] and Klar m Since then there have been many new developments in the construction 
of AP schemes for a large class of kinetic equations (cf. reviews by Degond [5] and Jin H3I). Time 
discretizations usually need an implicit-explicit (IMEX) approach [H [TUl3U I], exponential integration 
methods [7j, BGK-type penalty methods |8j, or micro-macro decomposition-based schemes (25[[27]. See 
also [28j l3i) . One key idea of the schemes is to split the equation into a nonstiff part, which is treated 
explicitly, and a stiff part, which will be implicit but can be implemented explicitly. The splitting 
should be taken in a way such that the combination preserves the AP property. Another possibility 
to treat the time-dependent case, which has been used often in the neutron transport literature, is 
to use an implicit semidiscretization in time first (e.g., backward Euler). This reduces the problem 
to the successive solution of steady problems (with a modified absorption term). For these, one can 
use a variety of techniques, among them the second order self-adjoint form of the transport equation. 
For a recent example, see M , where different spatial discretizations are investigated. The second 
order self-adjoint form, however, does not exist for the fully time-dependent, undiscretized transport 
equation. 

In this paper, we present a two-dimensional AP scheme for the time-dependent transport equa¬ 
tion 0 , where we do not use a semidiscretization in time. We combine a well-known scheme m which 
is based on the parity equations (which in turn are well-known [26]) with staggered grids, which in 
one dimension are fully understood m- As we have mentioned above, the advantage of the staggered 
approach, compared with the regular grid approach in [19] . is that in the diffusion limit we approach 
a compact stencil, as pointed out by Jin [13]. To be more precise, in one space dimension, using a 
regular Cartesian grid, the discrete diffusion limit of the scheme of [19] approximates the diffusion 
operator in 0 by (pi +2 ~ 2 pi + pi- 2 )/( 2 Ax) 2 (for the case cq = 1), while the current scheme gives the 
compact discretization (pi+i — 2 pi + pi- i)/(Ax) 2 which offers a better spatial resolution. 
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However, it is not trivial to extend the scheme from m to staggered grids in two dimensions. 
The first novelty of this paper is that we propose a staggering in two dimensions that requires fewer 
unknowns than one could have naively expected. The scheme shares similarities with diamond dif¬ 
ferencing but turns out to be different. The second contribution of this paper is that we present a 
rigorous stability analysis of the scheme m in one dimension and show uniform stability. Previously, 
only a stability argument was available. Similar to E3 EZj, we obtain an explicit CFL condition, 
which couples a hyperbolic and a parabolic condition and guarantees uniform stability. The scheme 
contains a splitting parameter that distributes terms between the explicit and implicit parts. For the 
choice of this parameter, our stability analysis yields an upper bound that is more restrictive than the 
one that has previously been used m- We discuss the choices of the methods we combine in the text. 
The advantages and disadvantages of these have been described in the original papers. For the sake 
of completeness, we repeat them in the text. 

The remainder of the paper is organized as follows. In section[2j we first derive the parity equations 
for the linear transport equation. Then, we describe the numerical method. In section [3j we show 
the AP property. We consider the asymptotic limit of the scheme and we state and prove a stability 
result, which gives a CFL condition and an upper bound on the relaxation parameter. Finally, several 
numerical tests are presented in section [4] to confirm the AP property of the scheme. Before the 
conclusion, we comment on several possible extensions (section [5]) . 


2 The Numerical Method 


First, we reformulate the transport equation into a parity equation. Then, we describe the angular, 
the spatial, and the time discretization of the method in detail. 

We begin with the transport equation in the diffusive scaling ([Tj) , restricted to two spatial dimen¬ 
sions x = ( x,y ). Following [19], the equation is split into four parts according to the quadrants of the 
velocity space. We obtain four equations with non-negative £, 77 > 0. This system can be rewritten if 
we define the even and odd parities 

r (1) (£,7?) = + /(-£> V)], r( 2 )(£,?7) = |[/(£,r/) + /(-£,- 77 )], 

j (1) (^r?) = ^[/(^,-7?) -/(-^ry)], j (2) (£, rj) = ^[/(£, rf) - /(-£, -rj)], 


leading to 

d*r« + £cy (1) - h<9 y j (1) 

<9 t r (2) + £<9 x j (2) + 7?<9 y j (2) 

+ 4^r« - ^a„r« 

+ Id,rW + 


^-(r (1) - p) - cr a r (1) + Q, 

£ z 

(4a) 

(r (2) - p) - cr a r (2) + Q, 

£ z 

(4b) 

£ z 

(4c) 


(4d) 


where p = ^ Jj , =1 fdv. It may seem that we have quadrupled the number of unknowns. Note, 
however, that due to symmetry we need to solve these equations for £, 77 in the positive quadrant only. 
Thus the number of unknowns is effectively the same. 

The discretization of the angular variable uses Gaussian quadrature points. Let 


£(A) = cos(A7r/2) and 77 (A) = sin(A7r/2) 
for every 0 < A < 1. Then, the density 

P=\ f [r ( 1 ) (£,h) + i' (2) (?>^)]^ A 


(5) 


0 


( 6 ) 
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can be approximated by a Gaussian quadrature rule on [0,1], where the quadrature points {A*} are 
mapped to {£j} and { 77 *} by <©■ 


2.1 Spatial Discretization 

First, we describe the spatial discretization, while keeping the time t continuous. We define a standard 
regular mesh with mesh size Ax x Ay and place the unknowns in the following way (see Figure [Taj): 

• r ( 1 ) ; it 2 ), p , and Q are located at the vertices (i,j) and the cell centers (i + 5 , j + 5 ); 

• jd) and j( 2 ) are located at the face centers (* + j) and ( i,j + 7). 

This choice enables us to approximate all spatial derivatives in the system Q by half-grid centered 
finite differences and yields a closed system of equations. 


2/7+1# —♦— # —*- 

Vji > —♦— < > —♦— t 

yj-ii —♦— —«- 


(a) Staggered grids 




(c) control volumes of ♦ 


Figure 1: Staggered grids and control volumes: red circles (vertices and cell centers), r^ 1 ), r^ 2 ), p , Q ; 
blue diamonds (face centers), j) 1 ), j) 2 ). 


In detail, the semidiscretized equations are defined as follows. Let i,j 6 Z. Then, the parities r^ 1 ' 
and r^ 2 ) satisfy (4a) and (4b), which on the vertices (i,j) and the cell centers (i + \,j + - 5 ) are given 
by 


](!) _:(!) 

a (!) , A i+ \,j J +i+5 a s, (1) \ (1) , ^ 

+ Z -- V -At/- = - Pi ^ - > 


Ax 


a (1) 1 >"^+1,7+5 ‘**>7+5 '**+5,7+1 "**+§,j < 7 s / (1) \ 

+ «-fc-’ - ” Ay = ~ P ‘+yK 

1 + Qi+i.i+i > 


;(!) 


■jW 


J 


Ay 

i(!) 


j (1) 


-a r (1) 

2 ,J 1 2 


(7) 


and 


.( 2 ) _ ( 2 ) .( 2 ) _ ( 2 ) 
a ( 2 ) , J *-5,7 . J +7+| 

MJ + z Ax +v- 


— — — (y-t 2 ') — 0 . .)_ a r ( 2 )+n. . 

— 2 Pi ,]) aaL i,j ^ hh,7 ! 


a _( 2 ) 

1 i+y+h 


Ay 

:( 2 ) _ • ( 2 ) .( 2 ) _ -( 2 ) 

J i+l,j+| *,7+5 . ‘+5,7+1 ++5,7 0's/ (2) s 

f-1-- + T -- = -sL-A.i+l -Pi+Jj+i) 


Ax 


Ay 


e 2 v 

( 2 ) 

- <r ++i j+i + i+i 


(8) 
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Similarly, the equations for the parities j -^ and on the face centers (z + \,j) and ( i,j + are 
derived from the equations (4c) and (4d): 


' r 


(i) 


(i) 


4 r i+ l,j T i,j 1 _ <?t -(1) 


(1) 


— r 


(i) 


' e 2 Ax 

,(i) 


Ay 


,(i) 


o •(!) 4 i+jj+j V T i,j +1 r ZJ 

e 2 Ax e 2 Ay 


£ 2J i+| J ’ 


^:(1) 

£ 2j *.i+l 


(9) 


and 


aft, + 


,W 

r *+ij 


£ r„- „■ - r, 


( 2 ) 

hi 


„( 2 ) 


-( 2 ) 

1 • • 1 • 1 


' £ 2 Ax 


( 2 ) 


n i+jj+i = -(2) 

£ 2 Ay £ 2J *+|d ’ 

( 2 ) 


t- i , , — r' i -W _ A 1 ) 

a .(2) , 4 >+2,j + j *- 5 J+2 , V T i,i +1 r i,i _ Cr t;(2) 

C'ij • • i 1 H-9 -7-1-9 -7- —-9 J ■ • , 1 


hi + 2 


Ax 


Ay 


£ 2J m+! • 


( 10 ) 


Remark 1. The method can be interpreted as a finite volume method and is therefore conservative. To 
see this, note that all the unknowns are given on two regular grids: vertices and cell centers or vertical 
cell centers and horizontal cell centers. Merging the corresponding grids, we obtain two staggered 
nonregular grids as shown in Figure la. These can also be interpreted as control volumes of a finite 
volume method. The control volumes ()i,j are defined as shown in Figures lb. lc, where Oij is the 
diamond around the point ( Xi,yj ). 

We define the volume averages corresponding to the control volumes by, e.g., 


-W - = 

M ’ 


If 1 

/ r (i) d(x,y) with |0i,j| = -AxAy 
I vi ,j | JO; ,' ^ 


( 11 ) 


and integrate the system Q over the control volumes. For instance, (4a) integrated over the volume 
0i,j is given by 


d Al + 1 0 . 


J ^ (4<9zj (1) + hdy] {l] ) d{x, y ) = — (r|j - p iyj ) + Q itj 


The integral can be simplified using Gauss’s theorem 



(4<9 T j (i) +y<9yj (1) ) d(x,y) 




d(x,y ) 




d(x,y) , 


( 12 ) 


(13) 


where n is the outer normal vector. It remains to compute the integrals over the four edges of the 
diamond. These terms are approximated by the trapezoidal rule, e.g., the upper right part is then given 
by 


rW>y j+ i) 
'O i+ l >%') 


j (1) • n ) d(x,y) = 


< x ^v j+ 1) 


' (x. 

v r- 


u y i) 

^ 1 
' 2 




Ay 


Ax 2 +A y 2 
Ax 


d(x,y) 


yj Ax 2 +A y 2 


• (1) + 1 ( 1 ) 

J m+ 1 *+§,j 


. (14) 
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In total, we obtain the approximation for the integral 


AxA y 


'dOi 


j ( 1 ) -»)d(*,»)» A(jW j 


i (1) 

J *-§, j 


+ 


Ay 


i (1) , 


i (1) , 

^i,j— 2 


(15) 


9 tr (1) + — 
* Ax 


and therefore the same semidiscrete equation as before 0 

+ i~ y (C§ ’Cl) = -S (r S “ ^ + «W ' <“> 

In the same way, we obtain the semidiscretized equations ([8]) -(10). 

In the case of a Cartesian mesh Ax = Ay, the finite volumes are rotated by f5 degrees with respect 
to the axes, which seems to be odd. However, it seems that the choice of grid points and therefore the 
volumes are unique in the strategy we have adopted. In fact, it seems quite natural when we consider 
the fact that the quantity fl 2 ' 1 describes the number of particles in £,rj > 0 and £, rj < 0. Thus a flow 
into the diagonal direction makes sense, and therefore interfaces at f5 degrees to the axes. This will 
be investigated further in future work. 

Remark 2. Note that ([3]) can formally be inverted to obtain the density of the transport equation 

r(1) (|£|>M) +£ «0"(Oj (1) (l£l>M) f or £v < 0 » 
r(2) (l£l,M) + £ **0"(Oj (2) (lfl>M) f° r ^v > o • 




(17) 


However, in the numerical scheme the parities r* 1 ), r^ 2 ) and jW, are not given on the same 
spatial grid and need to be interpolated. This is similar to other approaches which are based on parity 
decompositions. 


2.2 Time Discretization 

For simplicity, we consider again semidiscretized equations. This time, we keep the spatial variables 
x and y continuous and apply the time discretization technique from m- The idea is to introduce a 
relaxation parameter f) = such that we obtain a linear hyperbolic system with stiff relaxation. 

Then the linear hyperbolic part, which is nonstiff, can be separated from the stiff relaxation step. 
First, we rewrite the system of © as the diffusive relaxation system 

<%r (1) +£<9xj (1) - r]dy] (l) = -^(r (1) - p) - cr a r (1) + Q , 

<%r (2) + ^j (2) + r]dyi {2) = (r (2) - p) - cr a r (2) + Q , 

£ 1 (18) 
+ 0£d x r (1) - cfrjdyx^ = [a s j (1) + (1 - £ 2 (f)fd x r (l) - (1 - £ 2 4>)yd y r (1) ] , 

<9ij (2) + <f€d x r (2) + t; brjdyr (2) = [cr s j (2) + (1 - £ 2 </>)£<9 x r (2) + (1 - er 2 0)?7«5» y r (2) ] 

with 0 < f A 1/e 2 - The condition (j) > 0 is necessary for the hyperbolicity of the left-hand side, 
whereas the condition <fi < 1/e 2 ensures that the bracketed term on the right-hand side has a well- 
defined limit for e —> 0. In Remark [3] below we comment on the role of f and how this choice differs 
from Klar’s scheme EH. 

Thus, we split the equation into two parts, the transport step, 

<9t ( 1) + ^x.j (1) - = -a a r« + Q , 

d t r {2) + fd x p + r,d y )W = - aa ^) + Q , 

<9tj (1) + </^ x r (1) - cfrjdyi^ = -<7 a j (1) , 

<9tj (2) + </£<9x-r (2) + (fr/dyi^ = -<r a j (2) , 


(19) 
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and the relaxation step, 

3 f r (1) = -^(r (1) - p) , 
e- 

d t r {2) = (r (2) - p) , 

= -^2 ksj (1) + (1 - e 2 (t>)£,d x r (1) - (1 - £ 2 </>)t?<V ( 1) ] , 
<9ti (2) = Kj (2) + (1 - £ 2 </>)£<V (2) + (1 - £ 2 </>)r?<9 y r (2) ] . 


Finally, we apply the explicit Euler method to the first step and the implicit Euler method to the 
second step. Note that the implicit Euler method can be implemented explicitly, since p is preserved 
in the second step (which can be seen by adding the first two equations). 

The fully discrete scheme is just splitting 0-|To| into the two steps (19)- (20). 


Remark 3. Klar [21 ] developed a similar decomposition, which corresponds to f = 0 in our framework. 
As mentioned in m, there are two major differences. First, the resulting system is only weakly 
hyperbolic and therefore well-posedness is an issue. Second, computations are still performed on the 
hole velocity space. Using the symmetries of the parities, the computations can be performed on the 
first quadrant and the computational cost can be reduced. 


Remark 4. To generalize the problem from a two-dimensional “flatland” to a full two-dimensional 
problem, the domain of the angular variable changes from the unit circle to the unit disc. In order 
to do this change, the Gauss integration rule needs to be substituted by a two-dimensional integration 
rule on the unit circle. 


3 The AP property 

In this section, we analyze the AP property of the above scheme in two steps. First, we derive the 
discrete asymptotic limit. Second, we analyze stability. 


3.1 The Diffusion limit 


In the same way as above, we consider the spatial and the time discretization separately. The limit, as 
e —> 0, of the time discretization is derived in [19j . Hence, it remains to investigate the discrete limit 
of the spatial discretization. To this end, we consider the diffusive limit e —> 0 of the semidiscretized 
equations 0-©. 

First, the limit of 0 and 0 for the parities r^ 1 ) and is given by 


r (1) - p . • 
l i, 3 ~ ri,j ) 

id+i = Pi +h,3+b ’ 


„( 1 ) 


i+n,3 + 


i (1) 


i (1) = 


.( 1 ) 
h+ij 


£ r„-,— r. 


(i) 




Ax 


+ 


r « 

V D\,j+\ 


— r 


(i) 


(J, 


Ay 


„(i) 


£ i +b,j+ 


<?t 


l — rf^i . t r ( 1 ) _ r ( 1 ) 

; *~2’J+2 _|_ ^ hj + 1 T j,j 

Ax (J s Ay 


( 21 ) 
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Inserting these equations into ([T]) yields 

£ Pi+lj" 2 Pij -|- Pi— ] j 


d t pu 


+ 


d tP i+ l 


<7 t (Ax) 2 

2£ri Pj+l,j+i ~ Pj+ij-i ~ Pi—|,j+| + 
at AxA y 

r Pi,j +1 — 2 pij + Pi,j-\ _ _ „ , r\ 

a t (Ay) 2 “ aapi ’ j + Vij ' ’ 

e 2 Pi+i 3 +\~ 2 Pi+\p+\ + pi-y+\ 


»+ 2 j‘+ 2 (Jf (Ax) 2 

2£ r 7 Pi+lj+1 — Pi+lj — Pi,j+1 + Pij 


+ 


at AxA y 

v 2 p i+\,i +f “ 2/ M>5 


O'* 


(Ay) 5 


= -o’, 


( 22 ) 


+ ®i+b,j+h 


Treating (|8j) and (10) in the same way as above, we additionally obtain the following differential 
equations for p: 


dtPij 


Pi+l,j 2 Pij + Pi-1,j 

a s (Ax) 2 

2 €pPi+l,j+l ~ ~ ^jj+i + 

a s AxA y 

p 2 Pi,j+l — 2pi,j + Pi,j -1 _ , 

— u i j 




*+ 1 J+2 


0's (Ay) 2 

2 ’ J 1 2 (J s (Ax) 2 

2£y Pi+lj'+l — Pi+l,j — Pij'+l + Pij 
CT. S AxA y 

P 2 p *+|d+l “ 2p *+!-5 + p *+ld-| 


= —c, 


o s (Ay) 2 

Adding up the equations, the middle terms cancel and we obtain 

ti 


Pi - ( + ^+SJ'+5 


9*tPi,j 


£ Pi+lj 2 pi j + pt—\ j y Pi,j -\-1 2pjj T pij—i 


Oc 


(Ax) 5 


o s 


(Ay) 5 


— OaPij T Q. 


*J > 


StPi 


*+ 2 >•?+ 2 (j s 


g 2 ^+§,J + | 2/3 i+|j + ^ + 

(Ax) 2 

V 2 p i +\, j+l ~ 2 p i +\,\ + Pi +\ j-k 


a. 


(Ay) 5 


= -o. 


^i+ii+1 +< 5i+|j+l 


(23) 


(24) 


Integrating over £ 2 + y 2 = 1 yields the semidiscretized diffusion equations on the vertices and the cell 
centers. Note that integrating (22) or (23) over £ 2 + ?y 2 = 1, the middle terms cancel as well and we 
get the same result. 

As expected, the spatial discretization with staggered grids leads to a compact five-point stencil 
for the diffusion equation ([2]) . Together with the results m on the limit of the time discretization ( |19| ) 
- (20), this also shows that the formal limit of our scheme coincides with the diffusion equation. 























ASYMPTOTIC PRESERVING TWO-DIMENSIONAL STAGGERED GRID METHOD 


9 


3.2 Stability 

We limit our discussion to the one-dimensional case (see Remark [9] for the two-dimensional case) and 
show uniform stability with e using the von Neumann analysis IMEZIES]- 

In the following, we consider the transport equation in slab geometry and assume that the cross 
section at = cr s + e 2 a a > 0 is independent of x E M (see Remark [ 8 ] for space-dependent scattering). 
Further, we consider a source-free two velocity model v E {±1}- Then, the even and odd parities 

= \[f{t,x,l) + f(t,x,-l)] and j(i, x) = x, 1) - /(*, x, -1)] (25) 


fulfill 


d t x + d x ] = -a a i , 
dt] + hd x r = -4cr s j - cr a j 


(26) 


and the numerical scheme has the following update rule: For k = 0,1, 2,. 


r k+ 2 = r fc - At{D x f + a a v k ) , 
jfc+| = jfe _ At(cf)D x r k + a a i k ) , 


r fc+1 = r fc+ 2 


•k+l 


(27) 


J = 


e 2 +cr. 


,AtJ 




2 — 


At 


e 2 +a s Af 


(1 - e 2 4>)D x t 


k+l 


where D x denotes the half-grid centered finite difference of the spatial derivative. We place r on the 
half grid points (m + ^)Ax and j on the full grid points mAx. For this scheme, we do not expect 
positivity, but we seek a uniform CFL condition. 

For a von Neumann analysis of the scheme, we expand the parities in Fourier series: 


OO OO 

r (x,t) = ^ ai(t)e lix and j (x,t) = ^ bi{t)e^ x . 

£=—oo £—OO 


(28) 


As no mixing between the Fourier modes occurs during the update of the solution, it is sufficient to 
consider the evolution of 


r (x,t) = ai(t)e^ x and i(x,t) = bi(t)e 


iix 


(29) 


for some I and to determine the growth factor matrix of the Fourier coefficients. First, we note that 
the staggered grid derivatives can be rewritten as 


(D x r) (h (m + 3 ) , t) = apt) 


0 ilh(m+ 1 ) giihm 

h 


= 2 |sin {f)e m{m+ ^ai(t) , 


ith(m +_ iCh(m .— A) 

(D x j)(hm,t) = b t {t )- - -= 21 sin (f ) e Uhm b e (t) , 


(30) 


with h := Ax. To shorten the notation, we define di := ^ sin (^). Then, the first update step of the 
Fourier coefficients is given by 


^k 

bk 


{t + At) 


1 — a a At —A tdi 

—A t(f>di 1 — a a At 



^k 


pk_ 


it) 


( 31 ) 
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and the second step is given by 


a-k 

h 


(t + At) 


At 

e 2 +u s At 


1 

(1 - £ 2 <f)di 


e 2 +<r s At 





bk_ 


(t) 


= :G 2 


(32) 


Thus, the growth factor matrix is G := G 2 • G± 


G = 


1 — a a At —A tde 

- e 4*£a l K A t(l - e 2 4>) + 1) £ v + ^ Af (g 2 ( 1 - ct a At) + dj At 2 (l - e 2 ^)) 


For stability, the eigenvalues of the matrix G are of main interest. They can be written as 

Ai, 2 = g ± \Jg 2 ~ det(G) 


(33) 


(34) 


with g being the half trace and det(G) being the determinant of G: 


1 1 


j 2 £ 2 +<X s A t 


(Af 2 d 2 (l — £ 2 <t>) + (1 — a a At)(2e 2 + o s At)) and 


det(G) = s 2 +a s A^ (( 1 - a aAt) 2 - cj>d 2 At 2 ) . 

Proposition 1. Let the time step At and the relaxation parameter 4> satisfy 

At < min |^-,ma x{^£h, \h 2 o t }| 


(35) 


and 


, f hot < 2e , 

0 < f < < Y 

I ^ 2 , otherwise. 


(36) 


(37) 


Then, the numerical scheme is L 2 -stable. 


Remark 5. Note that the three terms in the time step restriction (36) can be interpreted separately. 


First, the ^eh term comes from the advection operator. Second, the \h 2 ot term corresponds to the 
Courant limit for explicit diffusion. Third, the — term is a result of the explicit treatment of the 
corresponding relaxation term. In this case, we assume that this term is small, so that At < A is not 
the restrictive term. Otherwise the term could be treated implicitly, which would remove the restriction. 


Remark 6. Note that the above restriction on 4> is stricter than the one suggested in m 0 < 
(f < 4j. Moreover, the condition hat < 2g is satisfied, if and only if the hyperbolic condition 
At < max{|e/i, \h 2 at} = \£h holds. This means, there are the following two cases: 

e/i j 

h 2 a t 

In addition, as £ —)• 0 the time step restriction becomes At < min{4^, ^h 2 at}, which does not vanish. 

In the proof of the proposition, we use the von Neumann analysis. A complete overview of these 
stability conditions can be found in the lecture notes by Trefethen [33] . 


and 0 < f , 

and 0 < (f < \ . 


(38) 


hat A 2g : At < min 
hat > 2e : At < min 
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Proof. Stability follows from the von Neumann condition if we can show |Ai ; 2 | < 1 for Ai 7 ^ A 2 and 
|Ai 5 2 I < 1 for Ai = A 2 . To show these inequalities, we consider three different cases: two complex 
eigenvalues, two real eigenvalues, and one eigenvalue. Since g and det(G') are real-valued (35), the 
cases are equivalent to: g 2 < det(G), g 2 > det(G'), and g 2 = det(G'). 

Case g 2 < det(G) (two complex eigenvalues): If the eigenvalues Ai^ are complex, their real part 
is g and their imaginary part is ±iy/det(G) — g 2 . Thus, the stability condition |Ai t 2 1 2 < 1 is satisfied 
if det(G) < 1. For the determinant we have the following estimate 


det ( G ) = sMasAt ^ 1 ~ - 4>djAt 2 ) < £ >^ sAt (l - (Ta^t + 


(39) 


where we used that — d 2 = sin 2 (^) < A and the CFL condition At < d-. It remains to show that 

the last term of (39) is bounded by 1. This is equivalent to 


<<?s+D 


&a — 1 


(40) 


which in turn is satisfied under the condition At < maxjieti, \h?at} and the assumption (37). This 


is one of the reasons for the choice of the upper bound of 4> in the assumption (37). 

Case g 2 > det(G) (two real eigenvalues): The determinant of G is always positive and therefore 
the eigenvalues are either both positive or both negative, and their sign changes with the sign of g. 
Thus, it is sufficient to show Ai < 1 if g > 0 and A 2 > — 1 if g < 0. In particular, one can show that 
this is equivalent to 

det(G) + 1 + 25 > 0 . (41) 


The first inequality is generic 

det(G) + 1 - 2 g = ~ t {a 2 a e 2 + a s a a - d e ) > 0 , 


(42) 


since d 2 < 0. Whereas, the second inequality requires the CFL condition (36). More precisely, under 
the condition 0 < 4>e 2 < 1 and At < 


1 

cr a ' 


we obtain 


det(G) p\p2g>l + 2g— 1 + ^ (Af 2 d 2 (l — e 2 cf) + (1 — cr a At)(2e“ + a s At)^j 

— e^+lsAt( £ ~ ~ AW + a sAt) . (43) 

On the one hand, this is obviously nonnegative under the condition At < le/i. On the other hand, 
the second term can be rewritten as 

e 2 _ 4 A+ + asAt = e 2^ _ aa At) + At(a t - jjr) , (44) 

which is nonnegative under the condition At < A- and At < ^h 2 at. Together, this yields the desired 
inequality det(G) + 1 + 2g > 0. 

Case g 2 = det(G) (one eigenvalue): The eigenvalue of G is Ai = A 2 = g. Thus, we need to 

show \g\ < 1. But as det(G) + 1 + 2g > 0 and det(G) < 1 (see above cases) already imply g > —1, it 
remains to show g < 1. Since at = a s + e 2 a a > 0, at least one of the terms a a At, a s At is positive and 
we obtain 

At)(2£ 2 + a s At)) 

11 / o \ (4o) 

<c \ £ 2 +asAt ((! _ a aAt + a a At){2e 2 + a s At + a s At)) = 1 . 

□ 
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Remark 7. If there is neither scattering nor absorption and the conditions (36) and (37) hold, then 


the relaxation parameter satisfies <f> = 0 and the hyperbolic condition is always satisfied. Further, the 


A t 2 d? 


determinant, det(G') = 1, and the half trace, g = 14 —, coincide only if dj = — ^sin 2 (y) = 0, 
which is equivalent to ^ € 7rZ. In most cases, this does not occur and therefore the case g 2 = det(G) 
does not arise. Then, we obtain g 2 < det(G) and the eigenvalues are distinct and satisfy (Ai^l = 1, 
so that stability follows. 

Remark 8. If the cross sections are space-dependent, the above analysis is not valid. In practice, the 
CFL condition is replaced by a worst-case condition. This means that we replace a a and at in (36) 


and (37) by its maximum and minimum, 


CTa.max = max<T a (x) and a t . m \ n = min a fix) , 


(46) 


respectively. 

Remark 9. In two dimensions, one can expect that the stability result from Proposition^ carries over 
with the following changes. We replace h = min(Ax, Ay) and add a factor of \ in front of the time 
step to account for the presence of growth rates in each of the two spatial dimensions. 


4 Numerical Results 


In this section, we consider different numerical test cases to demonstrate the performance of our 
scheme. Since we did not examine boundary conditions, we only consider examples, where the solution 
is compactly supported away from the boundary. We implemented periodic boundary conditions, so 
that there is no influence of any discretization of boundary values. 

The numerical calculations are performed using the two-dimensional scheme described in Section [2] 
with the stability conditions from Section 3.2. This means, we first choose the number of grid points 


(N x N) for the staggered grids corresponding to the test case. Then, we determine the maximal time 
step (cf. Proposition [T] Remark [8j and Remark [9]) 


At := 0.9 • \ 


mm 


{^T7Z’ maX ^ e/l ’ 3^*,min}} 


(47) 


and define the relaxation parameter 



hat < 2e , 
otherwise. 


(48) 


with h := jj, ay, m a Y := max x <7 a (x), and a t , m ; n := min x ay(x). The angular discretization uses a 
Gaussian quadrature with 16 points on the interval [0,1] for A. As the quadrature points are mapped 
to the directions ^ and g with ([5]), we obtain 16 points per quadrant. In all test cases, we compare the 
numerical solution on a grid where the parameter e is resolved to a grid on which it is underresolved, 
thus demonstrating the AP property. 

In the remainder of this section, we describe the test cases and the numerical results in detail. 
We consider four test cases to show different aspects of the AP property. First, we focus on the 
e-dependence and investigate the convergence order in different regimes. In the second and third test 
cases, there are large spatial differences in the cross sections. The second test case is continuous and 
rotationally invariant, whereas in the third test case the material cross sections and the source term 
are discontinuous. These two test cases intend to demonstrate the performance in multiscale problems. 
The last test case investigates the stability of the scheme-dependent on the choice of the relaxation 
parameter q A 
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4.1 Convergence order 

We examine the order of convergence with respect to the spatial variable. We expect first or second 
order convergence depending on the used CFL condition. If a hyperbolic condition is used, the time 
step is proportional to h. As the explicit Euler method is used for the time discretization, we cannot 
expect more than first order convergence in h. Whereas if the parabolic condition is used, the time step 
is proportional to h 2 . Then, the explicit Euler method predicts 0(h 2 ) convergence. Moreover, centered 
differences, which are used for the spatial discretization, are as well a second order approximation in 
h. Thus, we expect that the error is proportional to 0(h) when the hyperbolic condition is used, 
and 0(h 2 ), respectively, when the parabolic condition is used. To estimate the convergence order, we 
compute the £ 2 -error E(N) between the solution computed on a N x N grid and a reference solution. 
Using two different values N\ and IV 2 , we then estimate the convergence order by 

f n 2 = _ log(E(jV 1 ))-log(^(iV 2 )) , , 

Nl log(JVi) - log(AT 2 ) 


4.1.1 Method of manufactured solutions 

For the method of manufactured solutions (MMS), we first choose some function f{t,x,y,£,ri) and 
compute a corresponding source term and an initial condition, so that the chosen function is a solution 
of the transport equation. Let 


f{t, x, y, £, rj) = exp(-t) sin(2vrx) 2 sin(2vry) 2 (l + rj 2 ) 


(50) 


with (x,y) E [0, l] 2 . Further, let the scattering cross sections be given by u a = 0 and a s = 1. Then, 
the corresponding source term is given by 


Q(t> x, y, £, rj) 


d t f + ev ■ V x f 



(51) 


and the initial condition is given by 


f(t = 0, x, y , £, rj ) = sin(27rx) 2 sin(27ry) 2 (l + y) 2 . 


(52) 


We use the source term and the initial condition to compute a solution with the above scheme. For 
different grid sizes and different values of e, we compare the computed densities with the analytic 
density 

p{t,x,y ) = \ / fdv' (53) 

^ Jn 

at time t = 0.1. The results are shown in Figure [2] and Table [l] They confirm second order convergence 
in the parabolic case. In the hyperbolic case, the convergence order is even slightly higher than 
expected. 
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Figure 2: Convergence order (MMS): t 2 -error as a func¬ 
tion of the spatial resolution. Hyperbolic (filled markers) 
or parabolic (empty markers) CFL condition. 
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^16 

7^64 

-^32 

771128 

-^64 

771256 

-^128 

£0 

1.60 

1.50 

1.35 

1.23 

£1 

1.98 

1.93 

1.86 

1.76 

£2 

2.02 

2.00 ' 

1.99 

1.97 

£3 

2.02 

2.01 

2.00 

2.00 


Table 1: Convergence order (MMS): The 
term Ej^J is the convergence rate when go¬ 
ing from N\ x N\ to N 2 X N 2 grid points 
for a fixed mean free path = 10 _fc , 
k = 0,1, 2, 3. The dashed line indicates the 
switch from the hyperbolic to the parabolic 
condition. 


4.1.2 Gauss test 


We consider an example case with a smooth initial condition and isotropic scattering 

f(t = 0 ,x,y,v) = 4^10-2 exp(— f^-) for (x, y) € [—1,1] x [—1,1] , 
Q = 0 , (j t = o s = \ , a a = 0 , and e = 1, 1CT 1 , 10 -2 . 


(54) 


Then, we compute the density p at time t = 0.1 for different grid sizes and different values of e, so that 
the CFL condition (47) changes form hyperbolic to parabolic. As a reference solution, we use a highly 
resolved solution with 512 x 512 grid points. Tableland Figure [3] agree with the above assertion, 
showing first order convergence when the hyperbolic condition holds and second order, respectively, 
when the parabolic condition holds. 



Figure 3: Convergence order (Gauss test): f 2 -error as 
a function of the spatial resolution. Hyperbolic (filled 
markers) or parabolic (empty markers) CFL condition. 
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1.36 
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1.96 

1.37 

1.15 
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£2 

2.03 

2.01 

2.09 

2.41 


Table 2: Convergence order (Gauss test): 
The term is the convergence rate when 
going from N\ x N\ to IV 2 X IV 2 grid points 
for a fixed mean free path = 10 -fc , k = 
0,1,2. The dashed line indicates the switch 
from the hyperbolic to the parabolic condi¬ 
tion. 


4.2 Variable scattering 

In this test case, we examine the performance of the scheme, when the scattering is space-dependent. 
Compared to the previous test case, we fix the scaling parameter s and modify the scattering cross 
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section. Let 


f(t = 0 ,x,y,v) = exp(-f^) for 

£ = too ’ ( 5 = 0 > cra = 0, and 


/ x / \ fc 4 (c + ^(c-V2) 2 

ct{x,y) = a s {x,y ) = 


(z,y) € [-1,1] x [-1,1] , 

c = \Jx 2 + y 2 < 1 , 
otherwise. 


(55) 


Note that the total cross section a t (x,y) can be periodically extended to a C' 2 -function and at ^ ,y > 
ranges from 0 to 100. This wide range compared to the size of the domain causes strong variations of 
the solution, which are a challenge for numerical schemes. 

We compute the solution up to time t = e on two different grids. One of the grids underresolves 
the length scale e = ^gg (32 x 32 grid points) and the other one resolves it (512 X 512 grid points). 
Comparing the solution at different times (t = e; see Figure |4]) , we observe that the density, 

computed on the underresolved grid, matches the behavior of the density, computed on the resolved 
grid. 


4.3 Two material test 

The two material test case is a slight modification of the lattice test, which was proposed in [3|. It 
models a domain with different materials by discontinuous material cross sections and a discontinuous 
source term in space. 

In this problem, the computational domain is a 5 x 5 square. Most of the domain is purely 
scattering, except for some purely absorbing squares of size 0.5, which are distributed around an 
isotropic source in the middle of the domain 


Q(x,y) = 


1, (x,y)e[ 2,3] 2 

0 otherwise. 


(56) 


In the absorbing spots (cf. Figure 5a), the absorption coefficient jumps from 0 to 100, while the 
scattering coefficient jumps from 1 to 0. Thus, there are diffusive and kinetic regimes, although the 
scaling parameter satisfies e = 1. We obtain a rapid change of the solution at the transition zones, 
which may cause difficulties in the numerics. 

We compute the density up to time t = 1.7 on a coarse grid (64 x 64) and on a fine grid (512 x 512). 
The solutions are shown in Figure [5] Again, we observe that the solution on the underresolved grid 
resembles the solution on a grid that is resolved. In the case of the resolved solution, the oscillations 
near the beam edges are due to the angular discretization. They are the well-known ray effects for 
finite discrete velocity models (cf. [2] and references therein, as well as [241129)). 


4.4 Relaxation parameters and stability 

In our final test, we consider different relaxation parameters. Proposition [I] suggests an upper bound 
on the relaxation parameter (j), which in the hyperbolic case is more restrictive than in the parabolic 
case. We expect that in certain example cases our scheme becomes unstable if <j) is too large. 

Similar to the Gauss test, let 

f{t = 0 ,x,y,v) = 47r . 5 ^ 10 - a ex P(~ 4 . 5 x 10 -^ ) for (*»y) G [-1,1] x [- 1 , 1 ] , 
e = l, Q = 0 , at = a s = 1 , a a = 0 , N = 300 , and t = 0.36 . 


(57) 
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4 

3 


y 





Figure 4: Variable scattering: Density at | = 0.1 (first row), | = 0.5 (second row), and ^ = 1.0 (third 
row), computed on a 32 x 32 grid (first colum) ora512x512 grid (second column). The third column 
shows the density on a cut along y = 0. 



Figure 5: Two material test: (a) Geometry - source (orange), purely scattering at = a s = 1 (white 
and orange), purely absorbing a t = o a = 100 (black), (b) and (c) Density p at t = 1.7, computed on 
a 64 x 64 grid (b) or 512 x 512 grid (c). Logarithmic scaling, values are limited to seven orders of 
magnitude. 
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CL 



(a) </> = 7^ (b) (j> = 4r 

Figure 6: Stability: Density on a cut along y = 0, computed on a 300 x 300 grid up to time t = 0.36 
using different values of (f>. 


Then, we compute the density on a N x N = 300 x 300 grid up to time t 
relaxation parameters 


<^i — ffr — x 10 4 and (f >2 — p- — 1. 


0.36 using different 


(58) 


In the first case, the relaxation parameter satisfies the assumption of Proposition [l] and the solution 
is stable (see Figure [6a|), whereas in the second case, the assumption is violated and the solution 
starts to blow up (see Figure [6b]). As a consequence, the upper bound on the relaxation parameter in 
Proposition [I] can in general not be substituted by the less restrictive upper bound (j) < Jj- 


5 Extensions 

5.1 Boundary conditions 

In the numerical examples, we consider only periodic boundary conditions. However, one can use 
different approaches to implement, for instance, inflow boundary conditions. In the following, we 
discuss how to adapt two different methods mm to our scheme. Let the computational domain be 
given by X = [0, L x \ x [0, L y \ for some L x , L y > 0. Then, inflow boundary conditions on a continuous 
level are given by 

f(t, x, v) = f m (t, x, v) for n ■ x > 0, x £ <9X, (59) 

where n is the outer normal vector and / m describes the inflow. 

Similar to EUdniEoj, the solution of a kinetic half-space problem, describing the zeroth order 
diffusion boundary conditions, can be approximated to determine the outflow function f out (t, x, v): 

f(t, x, v) = / out (t, x, v) for n ■ x < 0, x € dX. (60) 

This determines completely the density / on the boundary and with this also the parities on the 
boundary. Considering, for instance, the left boundary x = 0, the boundary conditions of the parities 
r( 4 ), j - 1 - 1 are given by 

r (1) (£, rf) + ej (1) (£, r/) = /(£, -r?) = /“(£, -rj) , 
r (1) (£c?) ~£j (1) (f,f/) = = f° ut (-£,v) ■ 


( 61 ) 
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In one dimension, the parity r^ can be placed on the grid iAx, i = 0,..., N x with N x = kjj- and 


Ax 

jW on the staggered grid (i + ^)Ax, i = —1,..., N x . Then i ,( A lies on the boundary and j^ 1 -* can be 
interpolated using the ghost points — \ and N x + \ [2T] . This can be extended to two dimensions if 
we use a_ one-dimensional interpolation orthogonal to the boundary and add suitable ghost points (see 
Figure 


7). Then, the boundary condition (61) becomes for r^ 1 ) and j = 1,.. 


, N y — 1 and N y = 


r S + 1 1 j i 


•(i) 


. +j?>. 

2-? 2« 


= ffij and rfJ - ? 


d! 

2’ J 


+j ( , 1> 
2* 


= /o° 


out 

J 


and for j^ 1 ) and j = + \,... ,N — \ 


K r -lj +r S 


I pd 1 ) _ fin 
+ £ \j JOj 


and - 


,(i) 


T' 1 ■ 

2 ’-? 


_ rout 

£ h,j Jo,j 


(62) 


(63) 


Note that the edges of the boundary can be dropped, because these points are needed neither for 
computing the inner points nor for the interpolation (due to the one-dimensional interpolation rule). 
This choice yields a system with the same number of equations and unknowns. 

i i Q _ ~0 - i i 

L v i—-T- OOO —r—i 

-9 

o <>•♦•<> o 

Of-- J — 0 6 0 —~ ~ 1 
"---^--o--0 

U ljrp 


Figure 7: Black solid line denotes the boundary of the domain <9X; filled markers denote inner points; 
empty markers denote boundary and ghost points. 


Alternatively, the boundary conditions presented in m can be used. The key idea is to approxi¬ 
mate the j-unknowns on the boundary by 

jO = _ r/dyi 0) and = —^^r^ 2 ) — rjd y x ® (64) 

and insert this into the inflow boundary parts. For example, at the left boundary x = 0, the boundary 
condition for r^ 1 ! is given by 

/ m (£, ~v) = !' (1) (£, V) + £ J (1) (£, V) = r(1) (£, V)~£ (^r (1) + rjdyi^ . (65) 

To implement the boundary conditions in our scheme the spatial derivatives can, for instance, be 
approximated by one-sided finite-differences. There are multiple choices to do this, which need to be 
tested. 


5.2 Anisotropic scattering 

The simplicity of the scheme relies on the simplicity of the scattering operator, which is isotropic. In 
the following, we briefly consider the case of anisotropic scattering. Let the scattering operator be 
given by 


Kf(v) 


s(v , v')f(v')dv' 


( 66 ) 
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for some symmetric scattering kernel s, which depends on the cosine of the scattering angle v ■ v'. 

To compute the parity equations, the transport equation is rewritten for v = (£, 77), (£, —rj), (—£, r /), 
and (—£, —77) with nonnegative £, 77 > 0 and added in correspondence with the definition of the 
parities Q. This leads to a summation of the scattering operator of the form Kf(v) PKf(-v), which 
has to be expressed in terms of the parities. Since the scattering kernel satisfies s(—v,v') = s(v, —v'), 
we obtain 


K f(v) ± Kf(—v) = j s(v, v')f(v')dv' ± J s(—v,v')f(v')dv' 

= j s{v,v')f(v')dv' ± J s{y,v')f{-v')dv' = J s(y, v'){f{v') ± f(-v'))dv' . (67) 
More precisely, we consider for £, 77 > 0 

\{Kf{i,-rt) + Kf{-^ri)) = J s(£,ri,g,r/)^ (/(£', ~rf) + /(-£', rf)) d(g,rf) . (68) 

In this term, we recover the parities r^ and r^ 2 ) 


;(/(£', ~v') + f(-e,v')) = 


r^ for rf > 0 or £', 77' < 0 , 
r^ 2 ) otherwise. 


(69) 


Splitting the integral into the four quadrants and changing to positive 77', we obtain 


7 jKf(Z,~v) + Kf{-i,ri) = j ,rf) + a{£,ri,-g ,rf)d(g ,-rf) 

+ / + s(Z,ri,-g,rf))iW(g,r/)d(£ , ,T/) . (70) 

J £,i]>0 

Similarly, the other summations can be expressed in the parities. To simplify the notation, we define 
the following operators: 


V) : = / («(£, V, rf) + s(& V, -*/))/(£'> rf)d (£', rf) , 

•le,7?>o 

*/) := [ (s(£» rj, - 7 /) + s(£, 77 , rf))f{g, rf)d{tf , ?/) . 

Then, the parity equations for anisotropic scattering are given by 

d t r« + ^j (1) - 77^yj (1) = - J( r(1) - 77 + r (1) - Ifr^) - a a rW + Q , 

9 t r( 2 ) + £0 s j( 2 > + rjdy^) = - J(r( 2 ) - K+i<® - K~i<V) - <r a + Q , 

9d (1 > + W" - 4<V (1) = - K- jM) , 

d t j (2) + j^,r( 2 ) + ^,r< 2) = -^|(j( 2 ) - K+ j( 2 ) - JT- jW) - a a j( 2 ) . 


(71) 


(72) 


One could now apply a similar splitting as above, but in contrast to the isotropic case the implemen¬ 
tation of the relaxation step is not straightforward. It is necessary to invert the integral operators at 
a potentially expensive computational cost (cf. [T 7 j for more details). 
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6 Conclusions 

In this paper, we have introduced a two-dimensional AP scheme for the linear transport equation. The 
linear transport equation has the diffusion equation as an analytic asymptotic limit. For AP schemes 
the discretization has to be chosen such that the analytic limit is preserved at a discrete level and 
the scheme is uniformly stable with respect to the mean free path. Here, we used a parity-based time 
discretization combined with a staggered-grid spatial discretization. 

We have shown that the spatial discretization has the desired AP property. In particular, due to 
the use of staggered grids, a compact five-point stencil can be achieved in the limiting discrete diffusion 
limit. Furthermore, the parity-based time discretization is suitable for the use of staggered grids, as 
the coupling between the even and odd parities reduces the number of the required unknowns. In 
addition, we have presented a rigorous stability analysis for the same scheme in one dimension. This 
provides a condition on the relaxation parameter and a CFL condition. Finally, we have performed 
several numerical tests for the two-dimensional scheme, which demonstrate the AP property. Since 
staggered grids can easily be extended to three dimensions, there is a straightforward generalization 
of our method to three spatial dimensions. Although we did not test the method, we expect that it 
has similar properties. 

In the future, it would be worthwhile to investigate the time discretization. Since our method uses 
a simple time-integration method (explicit Euler method), the convergence order is in general limited 
to one. To maintain a second order scheme, one could use some higher order IMEX time-integration 
method. Another possible scope of future work is to apply staggered grids in combination with a 
parity-based time discretization to other kinetic equations. 


References 

[1] S. Boscarino and G. Russo. On a class of uniformly accurate IMEX Runge-Kutta schemes and 
applications to hyperbolic systems with relaxation. SIAM Journal on Scientific Computing , 
31(3):1926-1945, 2009. 

[2] T. Brunner. Forms of approximate radiation transport. SAND2002-1778, Sandia National Lab¬ 
oratory, 2002. 

[3] T. Brunner and J. Holloway. Two-dimensional time dependent Riemann solvers for neutron 
transport. Journal of Computational Physics, 210(l):386-399, 2005. 

[4] R. E. Caflisch, S. Jin, and G. Russo. Uniformly accurate schemes for hyperbolic systems with 
relaxation. SIAM J. Numer. Anal., 34(1):246-281, 1997. 

[5] P. Degond. Asymptotic-preserving schemes for fluid models of plasmas. Panoramas et syntheses, 
39-40:1-90, 2013. 

[6] J. Deng. Implicit asymptotic preserving schemes for semiconductor Boltzmann equation in the 
diffusive regime. Int. J. Numer. Anal. Model., 11(1): 1—23, 2014. 

[7] G. Dimarco and L. Pareschi. Exponential Runge-Kutta methods for stiff kinetic equations. SIAM 
Journal on Numerical Analysis, 49(5):2057-2077, 2011. 

[8] F. Filbet and S. Jin. A class of asymptotic-preserving schemes for kinetic equations and related 
problems with stiff sources. Journal of Computational Physics, 229(20):7625-7648, 2010. 

[9] F. Golse, S. Jin, and C. Lever more. The convergence of numerical transfer schemes in diffusive 
regimes I: Discrete-ordinate method. SIAM journal on numerical analysis, 36(5):1333~1369, 1999. 



ASYMPTOTIC PRESERVING TWO-DIMENSIONAL STAGGERED GRID METHOD 


21 


[10] F. Golse and A. Klar. A numerical method for computing asymptotic states and outgoing distri¬ 
butions for kinetic linear half-space problems. Journal of statistical physics, 80(5-6):1033-1061, 
1995. 

[11] S. Jin. Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms. J. 
Comput. Phys., 122(1):51-67, 1995. 

[12] S. Jin. Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. SIAM 
Journal on Scientific Computing , 21(2):441-454, 1999. 

[13] S. Jin. Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: A 
review. In Lecture Notes for Summer School on Methods and Models of Kinetic Theory(M&MKT), 
Porto Ercole (Grosseto, Italy), June 2010 , pages 177-216. Riv. Mat. Univ. Parma, 2012. 

[14] S. Jin and C. D. Levermore. Fully discrete numerical transfer in diffusive regimes. Transport 
Theory Statist. Phys., 22(6):739-791, 1993. 

[15] S. Jin and C. D. Levermore. Numerical schemes for hyperbolic conservation laws with stiff 
relaxation terms. J. Comput. Phys., 126(2) :449-467, 1996. 

[16] S. Jin and D. Levermore. The discrete-ordinate method in diffusive regimes. Transport Theory 
Statist. Phys., 20(5-6):413-439, 1991. 

[17] S. Jin and L. Pareschi. Discretization of the multiscale semiconductor boltzmann equation by 
diffusive relaxation schemes. Journal of Computational Physics, 161:312-330, 2000. 

[18] S. Jin, L. Pareschi, and G. Toscani. Diffusive relaxation schemes for multiscale discrete-velocity 
kinetic equations. SIAM Journal on Numerical Analysis, 35(6):2405-2439, 1998. 

[19] S. Jin, L. Pareschi, and G. Toscani. Uniformly accurate diffusive relaxation schemes for multiscale 
transport equations. SIAM Journal on Numerical Analysis, 38(3):913-936, 2000. 

[20] A. Klar. Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semi¬ 
conductor equations. SIAM Journal on Scientific Computing, 19(6):2032-2050, 1998. 

[21] A. Klar. An asymptotic-induced scheme for nonstationary transport equations in the diffusive 
limit. SIAM journal on numerical analysis, 35(3):1073-1094, 1998. 

[22] E. Larsen and J. Morel. Asymptotic solutions of numerical transport problems in optically thick, 
diffusive regimes II. Journal of Computational Physics, 83(l):212-236, 1989. 

[23] E. Larsen, J. Morel, and W. Miller Jr. Asymptotic solutions of numerical transport problems in 
optically thick, diffusive regimes. Journal of Computational Physics, 69(2):283-324, 1987. 

[24] K. Lathrop. Ray effects in discrete ordinates equations. Nucl. Sci. Eng, 32(3):357-369, 1968. 

[25] M. Lemou and L. Mieussens. A new asymptotic preserving scheme based on micro-macro formu¬ 
lation for linear kinetic equations in the diffusion limit. SIAM Journal on Scientific Computing, 
31(l):334-368, 2008. 

[26] E. E. Lewis and W. F. Miller. Computational Methods of Neutron Transport. John Wiley and 
Sons, 1984. 

[27] J.-G. Liu and L. Mieussens. Analysis of an asymptotic preserving scheme for linear kinetic 
equations in the diffusion limit. SIAM Journal on Numerical Analysis, 48(4):1474-1491, 2010. 



22 


KERSTIN KUPPER, MARTIN FRANK, AND SHI JIN 


[28] L. Mieussens. On the asymptotic preserving property of the unified gas kinetic scheme for the 
diffusion limit of linear kinetic models. J. Comput. Phys., 253:138-156, 2013. 

[29] W. Miller and W. Reed. Ray-effect mitigation methods for two-dimensional neutron transport 
theory. Nucl. Sci. Eng , 62:391-411, 1977. 

[30] J. E. Morel, B. T. Adams, T. Noh, J. M. McGhee, T. M. Evans, and T. J. Urbatsch. Spatial 
discretizations for self-adjoint forms of the radiative transfer equations. Journal of Computational 
Physics , 214(1) :12-40, 2006. 

[31] L. Pareschi and G. Russo. Implicit-explicit Runge-Kutta schemes for stiff systems of differential 
equations. Recent trends in numerical analysis, 3:269-289, 2000. 

[32] B. Seibold and M. Frank. StaRMAP-A second order staggered grid method for spherical har¬ 
monics moment equations of radiative transfer. ACM Transactions on Mathematical Software, 
41, 2014. 

[33] L. Trefethen. Finite difference and spectral methods for ordinary and partial differential equations, 
unpublished text, 1996. available at http://people.maths.ox.ac.uk/trefethen/pdetext.html. 

[34] K. Xu and J.-C. Huang. A unified gas-kinetic scheme for continuum and rarefied flows. J. Comput. 
Phys., 229(20):7747-7764, 2010. 



